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Abstract. The capacity to model magnetohydrodynamical (MHD) flows in 
dynamical, strongly curved spacctimcs significantly extends the reach of numerical 
relativity in addressing many problems at the forefront of theoretical astrophysics. 
We have developed and tested an evolution code for the coupled Einstcin- 
Maxwell-MHD equations which combines a BSSN solver with a high resolution 
shock capturing scheme. As one application, we evolve magnetized, differentially 
rotating neutron stars under the influence of a small seed magnetic field. Of 
particular significance is the behavior found for hypermassive neutron stars 
(HMNSs), which have rest masses greater the mass limit allowed by uniform 
rotation for a given equation of state. The remnant of a binary neutron star merger 
is likely to be a HMNS. We find that magnetic braking and the magnetorotational 
instability lead to the collapse of HMNSs and the formation of rotating black holes 
surrounded by massive, hot accretion tori and collimated magnetic field lines. 
Such tori radiate strongly in neutrinos, and the resulting neutrino-antineutrino 
annihilation (possibly in concert with energy extraction by MHD effects) could 
provide enough energy to power short-hard gamma-ray bursts. To explore the 
range of outcomes, we also evolve differentially rotating neutron stars with lower 
masses and angular momenta than the HMNS models. Instead of collapsing, the 
non-hypermassive models form nearly uniformly rotating central objects which, 
in cases with significant angular momentum, are surrounded by massive tori. 
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1. Introduction 

Differential rotation is likely to be the norm in many astrophysical settings. For 
example, differentially rotating neutron stars can form from the collapse of massive 
stellar cores, which acquire rapid differential rotation during collapse even if they 
are spinning uniformly at the outset ^ |5] (see also Differential rotation can 

also arise from the mergers of binary neutron stars |SJ In these new-born, 
dynamically stable, neutron stars, magnetic fields and/or viscosity will transport 
angular momentum and substantially change the configurations on a secular timescale. 

Some newly-formed differentially rotating neutron stars may be hypermassive. 
The mass limits for non-rotating stars [the Oppenheimer-Volkoff limit] and for rigidly 
rotating stars (the supramassive limit, which is only about 20% larger) can be 
significantly exceeded by the presence of differential rotation |H] . Mergers of binary 
neutron stars could lead to the formation of such hypermassive neutron stars (HMNSs) 
as remnants. The latest binary neutron star merger simulations in full general 
relativity OEIE] have confirmed that HMNS formation is indeed a possible outcome. 
HMNSs could also result from core collapse of massive stars. Differentially rotating 
stars tend to approach rigid rotation when acted upon by processes which transport 
angular momentum. HMNSs, however, cannot settle down to rigidly rotating neutron 
stars since their masses exceed the maximum allowed by rigid rotation. Thus, delayed 
collapse to a black hole and, possibly, mass loss may result. 

The merger of binary neutron stars has been proposed for many years as an 
explanation of short-hard GRBs |12II13| . According to this scenario, after the merger, 
a stellar-mass black hole is formed, surrounded by hot accretion torus containing ~ 1— 
10% of the total mass of the system. The GRB fireball is then powered by energy 
extracted from this system, either by MHD processes or emission from neutrino- 
antincutrino annihilation. The viability of this model depends on the presence of 
a significantly massive accretion disk after the collapse of the remnant core, which in 
turn depends on the mechanism driving the collapse. 

Though magnetic fields likely play a significant role in the evolution of HMNSs, 
the numerical tools needed to study this problem have not been available until 
recently. In particular, the evolution of magnetized HMNSs can only be determined by 
solving the coupled Einstein-Maxwell-MHD equations self-consistently in full general 
relativity. Recently, Duez et al. |2] and Shibata and Sekiguchi ^B] independently 
developed codes designed to do such calculations for the first time (see also 117) 1. 
The first simulations of magnetized hypermassive neutron star collapse (assuming 
both axial and equatorial symmetry) were reported in [T^J, and the implications of 
these results for short GRBs were presented in fP^. These simulations proved that 
the amplification of small seed magnetic fields by a combination of magnetic winding 
and the magnetorotational instability (MRI) (201 E] is sufficient to trigger collapse in 
hypermassive stars on the Alfven timescale, confirming earlier predictions [71 122|. 

We also compare the results for hypermassive stars with the evolution of two 
differentially rotating models below the supramassive limit in order to highlight the 
qualitatively different physical effects which arise in the evolution. Given a fixed 
equation of state (EOS), the sequence of uniformly rotating stars with a given rest mass 
has a maximum angular momentum J ma x- A nonhypermassive star having angular 
momentum J > J max is referred to as an "ultraspinning" star. We perform simulations 
on the MHD evolution of two nonhypermassive stars - one is ultraspinning and the 
other is not; we refer to the latter as "normal." Instead of collapsing, they evolve to 
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Table 1. Initial Models. 



Case 


M /M ,TOv a 


M /M , sup b 


R e JM c 


JIM' 1 d 


T mt /\W\ C 




A 


1.69 


1.46 


4.48 


1.0 


0.249 


0.33 


Bl 


0.99 


0.86 


8.12 


1.0 


0.181 


0.40 


B2 


0.98 


0.85 


4.84 


0.38 


0.040 


0.34 



a The ratio of the rest mass Mq to the TOV rest mass limit for the given EOS. 
b The ratio of the rest mass Mq to the rest mass limit for uniformly rotating stars of 
the given EOS (the supramassive limit). If this ratio is greater than unity, the star 
is hypermassive. 

c The equatorial coordinate radius i? cq normalized by the ADM mass. 
d The ratio of the angular momentum J to M 2 (the angular momentum parameter) . 
° The ratio of the rotational kinetic energy to the gravitational binding energy. 
f The ratio of the angular velocity at the equator to the central angular velocity. 

a new equilibrium state after several Alfven times. The normal star settles down to a 
uniformly rotating configuration. In contrast, the ultraspinning star settles down to a 
nearly uniformly rotating central core, surrounded by a differentially rotating torus. 

The key subtlety in all of these simulations is that the wavelengths of the MM 
modes must be well resolved on the computational grid. Since this wavelength is 
proportional to the magnetic field strength, it becomes very difficult to resolve for 
small seed fields. However, the simulations reported here succeed in resolving the 
MRI. In what follows, we assume geometrized units such that G = c = 1. 

2. Initial Models 

We consider three representative differentially rotating stars, which we call "A" , "Bl" , 
and "B2." Their properties are listed in Tabled Star A is hypermassive while stars Bl 
and B2 are not. These configurations are all dynamically stable. All three models are 
constructed using a T = 2 polytropic EOS, P = Kp^, where P, K, and po are the 
pressure, polytropic constant, and rest-mass density, respectively. The rest mass of 
star A exceeds the supramassive limit by 46%, while the rest masses of stars Bl and 
B2 are below the supramassive limit. The angular momentum of star Bl exceeds the 
maximum angular momentum (J max ) for a rigidly rotating star with the same rest 
mass and EOS, whereas star B2 has angular momentum J < J max - Thus, star Bl 
is "ultraspinning," while star B2 is "normal." These models may be scaled to any 
desired physical mass by adjusting the value of K 

Following previous papers (e.g, [7J |21 HH] ) j we choose the initial rotation law 
u°u v = A 2 (il c — ft), where it^ is the four-velocity, f2 c is the angular velocity along the 
rotational axis, and f2 = u^/u is the angular velocity. The constant A has units of 
length and determines the steepness of the differential rotation. For these models, A is 
set equal to the coordinate equatorial radius R e q- The corresponding values of f2 cq /i7 c 
are shown in Table (where fi oq is the angular velocity at the equatorial surface). 

We must also specify initial conditions for the magnetic field. We choose to 
add a weak poloidal magnetic field to the equilibrium model by introducing a vector 
potential of the following form A v — vD 2 ma.x[At,(P — P C ut),0], where the cutoff P cut 
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is 4% of the maximum pressure, and is a constant which determines the initial 
strength of the magnetic field. We characterize the strength of the initial magnetic 
field by C = max(S? x/(87rP)), i.e. the maximum value on the grid of the ratio of 
the magnetic energy density to the pressure (where refers to the magnetic field 
in the comoving frame of the fluid). We choose At, such that C ~ 10 -3 -10 -2 . We 
have verified that such small initial magnetic fields introduce negligible violations of 
the Hamiltonian and momentum constraints in the initial data. 

3. Numerical Methods 

Duez et al. 14 and Shibata and Sekiguchi |15| have independently developed 
new codes to evolve magnetized fluids in dynamical spacetimes by solving the 
Einstein-Maxwell-MHD system of equations self-consistently. Several tests have been 
performed with these codes, including MHD shocks, nonlinear MHD wave propagation, 
magnetized Bondi accretion, MHD waves induced by linear gravitational waves, and 
magnetized accretion onto a neutron star. Details of our techniques for evolving the 
Einstein-Maxwell-MHD system as well as tests can be found in |141 1TB) . We have 
performed several simulations for identical initial data using both codes and found 
that the results are essentially the same. The simulations presented here assume axial 
and equatorial symmetry. We adopt the Cartoon method for evolving the BSSN 
equations |27|. and use cylindrical coordinates for evolving the induction and MHD 
equations. In this scheme, the coordinate x is identified with the cylindrical radius zu, 
and the y-direction corresponds to the azimuthal direction. 

4. Results 

4.1. Star A 

We have performed simulations on star A with initial field strength given by C — 
2.5 x 10~ 3 . We use a uniform grid with size (N,N) in cylindrical coordinates (m,z), 
which covers the region [0, L] in each direction, where L = 4.5i? cq - For star A, 
i? cq = 4.5M = 18.6 km(M/2.8M0). To check the convergence of our numerical 
results, we perform simulations with four different grid resolutions: N — 250, 300, 400 
and 500. Unless otherwise stated, all results presented in the following subsections are 
from the simulation data with resolution N = 500. 

4- 1.1. General features of the evolution Figure ^ shows snapshots of the density 
contours and poloidal magnetic field lines (lines of constant A v ) in the meridional 
plane. In the early phase of the evolution, the frozen-in poloidal magnetic fields lines 
are wound up by the differentially rotating matter, creating a toroidal field which 
grows linearly in time (see Fig. EJi). When the magnetic field becomes sufficiently 
strong, magnetic stresses act back on the fluid, causing a redistribution of angular 
momentum. The core of the star contracts while the outer layers expand. As shown 
in Fig.0:, the effect of the MRI becomes evident for t ^ 6P C . This is manifested as a 
sudden increase in the maximum value of \B X \ (= |-B ro |), which grows exponentially 
for a short period (about one e-folding) before saturating. The effect of the MRI is 
also visible in Fig.^ where we see distortions in the poloidal field lines. The amplitude 
of the toroidal field begins to decrease after t J> 20P C ~ (see Fig. |2|I) and the core 
of the star becomes less differentially rotating. 
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Figure 1. Snapshots of rest-mass density contours and poloidal magnetic field 
lines for star A at selected times. The first and third rows show snapshots of the 
rest-mass density contours and velocity vectors on the meridional plane. The 
second and fourth rows show the corresponding field lines (lines of constant 
Atp) for the poloidal magnetic field at the same times. The density contours 
are drawn for po/Pmax(0) = lO" - 36 ^ 09 (i = 0-10), where p max (0) is the 
maximum rest-mass density at t = 0. The field lines are drawn for A v = 
+ (^,max - A iPtlnin )i/20 (i = 1-19), where A v ,max and A^, jIT1 i n are the 
maximum and minimum values of A v , respectively, at the given time. The thick 
solid curves denote the apparent horizon. In the last panel, the field lines are 
terminated inside the black hole at the excision boundary. 
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Figure 2. Evolution of the central rest-mass density p c , central lapse a c , and 
maximum values of \B X \ and \B V \. |B x | ma x and |-B y | ma x are plotted in units of 
\J Pmax(0). The solid, long-dashed, dot-dash, and dotted curves denote the results 
with 7V=250, 300, 400, and 500 respectively. The dashed line in (c) represents 
an approximate slope lj = 0.18/P C for the exponential growth rate of the MRI, 
5B X oc e ut . The dashed line in (d) represents the predicted growth of -B^lmax at 
early times from linear theory. 



The combined effects of magnetic braking and the MRI eventually trigger 
gravitational collapse to a black hole at t » 66P C ~ 36(Af/2.8Af Q ) ms, where P c 
is the initial central rotation period. A collimated magnetic field forms near the polar 
region at this time (see Fig. However, a substantial amount of toroidal field is 
still present. Without black hole excision, the simulation becomes inaccurate soon 
after the formation of the apparent horizon because of grid stretching. To follow the 
subsequent evolution, a simple excision technique is employed [2HI EH • We are able 
to track the evolution for another 300 M m 8P C . We find that not all of the matter 
promptly falls into the black hole. The system settles down to a quasiequilibrium 
state consisting of a black hole surrounded by a hot torus and a collimated magnetic 
field near the polar region (see the panels corresponding to time t = 74.6P C in Fig.^J. 

4-1-2. Resolution study Four simulations were performed with different resolutions 
(see Fig. EJ): N = 250, 300, 400 and 500. We find that the results converge 
approximately when N ^ 400. On the other hand, results are far from convergent 
for N SS 300. Most importantly, the effect of the MRI is not captured in the two 
lowest resolution runs, as shown by the behavior of l-B^lmax in Fig. [21 This is because 
the wavelength of the fastest growing MRI mode (A max ) is not well-resolved for low 
resolutions. We find that we need a resolution A/A max <, 0.14 (N ^> 400) in order 
to resolve the MRI modes. In contrast, the growth of B y by magnetic winding is 
resolved for all four resolutions. The straight dashed line in Fig. 01 corresponds to the 
prediction of linear theory for the growth rate due to winding. This slope agrees with 
the actual growth of |-B y | max in the early (magnetic winding) phase of the simulation, 
but as back-reaction (magnetic braking) becomes important, the toroidal field begins 
to saturate. (See jHJ f° r a discussion of the MRI and winding effects in linear theory.) 
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Figure 3. Evolution of the irreducible mass and the total rest-mass outside 
the apparent horizon. (Here, tah is the local coordinate radius of the apparent 
horizon.) 



4-1.3. Evolution with excision Soon after the formation of the apparent horizon, the 
simulation becomes inaccurate due to grid stretching and an excision technique is 
required to follow the subsequent evolution. For details on our excision techniques, 
see [201 CHH' During the excision evolution, we track the irreducible mass of 
the black hole by computing the area of the apparent horizon *4ah and using 
Mj rr ps yj _4ah/167t. The irreducible mass and the total rest mass outside the apparent 
horizon are shown in Fig. 03 The total ADM mass of the final state system, consisting 
of a BH surrounded by a massive accretion torus, is well defined. In contrast, there is 
no rigorous definition for the mass of the black hole itself. We obtain a rough estimate 
following the steps described in [31]. We find that M^oie ~ 0.9M, where M is the 
total ADM mass of the system, and Jhoie/^hoic ~ 0-8. 

The black hole grows at an initially rapid rate following its formation. However, 
the accretion rate Mq gradually decreases and the black hole settles down to a quasi- 
equilibrium state. By the end of the simulation, Mo has decreased to a steady rate 
of pa 0.01Mo/Pc, giving an accretion timescale of ~ 10-20P C pa 5-10 ms(M/2.8M ). 
Also, we find that the specific internal thermal energy in the torus near the surface 
is substantial because of shock heating. The possibility that this sort of system could 
give rise to a GRB is discussed in 

4-2. Star A, comparison of different magnetic field strengths 

In order to test the scaling of our results for different values of the initial magnetic 
field strength, we have examined three other values of C in addition to the value 
of 2.5 x 1CP 3 chosen for the results of Section 14.11 Namely, we consider C = 
{1.25,2.5,5.0,10} x 10~ 3 , and the results are shown in Fig. 01 For the portion of 
the simulations in which magnetic winding dominates, the behavior is expected to 
scale with the Alfven time [221- In other words, the same profiles should be seen for 
the same value of tjtj^. From the left panel of Fig.QJ it is evident that this scaling holds 
very well for the toroidal field and for the central density and lapse, while t < OAtA ■ 
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Figure 4. Left: Selected parameters plotted against scaled time (t/tA) for 
evolutions of star A with four different magnetic field strengths: C = 1.25 X 10~ 3 
(solid lines), C = 2.5 X 10 -3 (long-dashed lines), C = 5.0 X 10 -3 (short- 
dashed lines), and C = 10 -2 (dotted lines). All runs were performed with 
the same resolution (400 2 zones with outer boundaries at 20M). When plotted 
against scaled time, the curves line up at early times (t £ 0.54a = life) when 
the evolution is dominated by magnetic winding. Right: Maximum value of 
\B X \ plotted vs. t/P c for evolutions of star A with four different magnetic field 
strengths. The line styles are the same as in the left panel. The behavior of 
I B x | max is dominated by the effects of the MRI and thus does not scale with the 
Alfvcn time. The curves corresponding to the two highest values of C (dotted 
and dashed) terminate at the time when the star collapses 



The later evolution is driven mainly by the MRI, which does not scale with the Alfven 
time. The scaling also does not hold during the collapse phase, when the evolution 
is no longer quasi-stationary. Though the scaling breaks down at late times in these 
simulations, the qualitative outcome is the same in all cases. The behavior of l-B^lmax 
for these four different values of C is shown in the right panel of Fig. The sudden 
sharp rise of l-B^lmax signals the onset of the MRI, and the approximate agreement 
of the slopes for different values of C indicates that the exponential growth rate of 
the MRI does not depend on the initial magnetic field strength (as expected from the 
linear analysis). 



4.3. Star Bl 

Here, we present results for the evolution of star Bl with C = 2.5 x 1CP 3 . This run was 
performed with resolution 400 2 and outer boundaries located at 4.5i? eq (36. AM). Since 
this model is not hypermassive, the redistribution of angular momentum through MHD 
effects will not lead to collapse. However, since this star is ultraspinning and angular 
momentum is conserved in axisymmetric spacetimes, it cannot relax to a uniform 
rotation state everywhere unless a significant amount of angular momentum can be 
dumped to the magnetic field. We find that this model simply seeks out a magnetized 
equilibrium state which consists of a fairly uniformly rotating core surrounded by a 
differentially rotating torus. 

Figure presents the evolution of some relevant quantities for this case. From the 
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Figure 5. Evolution of central rest-mass density p c , central lapse a c , maximum 
values of \B X \ and \B»\ for star Bl. The magnetic fields l-B^lmax 

and \B»\ max are 

plotted in units of \J p m ax (0) . Note that the lines become fairly horizontal at late 
times, indicating that an equilibrium has been reached. The dashed line in (c) 
represents an approximate slope of u) = (0.37/P c ) for the exponential growth rate 
of the MRI, 5B X cx e ut . The dashed line in (d) represents the predicted growth 
of \B y \ max from linear theory. 



central density and lapse, it is evident that the star has settled into a more compact 
equilibrium configuration. This is consistent with the expectation that magnetic 
braking should transfer angular momentum from the core to the outer layers. A 
brief episode of poloidal magnetic field growth due to the MRI is indicated by the plot 
of | B x | max in Figure [SJ The instability saturates and quickly dies away, leaving the 
strength of the poloidal field largely unchanged. Early in the evolution, the maximum 
value of the toroidal component \B V \ rises due to magnetic winding. This growth 
saturates at t ~ 10P C ~ 0.5^ . We note, however, that the toroidal magnetic field is 
non-zero in the final equilibrium state, though it is no longer growing due to magnetic 
winding. 

Snapshots of the evolution in the x-z plane are shown in Fig. |5| The density 
contours for times t = through 25.0P C show that angular momentum redistribution 
leads to the formation of a more compact star surrounded by a torus. At t = 10P C , 
the distortions of the magnetic field lines due to the MRI are clearly visible. As 
the disk expands, magnetic field lines attached to this low density material open 
outward, eventually leading to the field structure seen in the last three times shown 
in Fig. El for which some field lines are still confined inside the star while others have 
become somewhat collimated along the z-axis. For t i> 35P C , the density contours 
and poloidal magnetic field lines change very little, indicating that the system has 
reached an equilibrium state which is quite different from the initial state. We note 
that a significant toroidal field persists at late times when the system has essentially 
settled down to a final state. In addition, the final configuration also has significant 
differential rotation, especially in the outer layers. However, this differential rotation 
no longer winds up the magnetic field lines (i.e., the toroidal field strength does not 
grow). This is because the rotation profile has adjusted so that Q is approximately 
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Figure 6. Snapshots of density contours and poloidal magnetic field lines for 
star Bl. The first and third rows show snapshots of the rest-mass density 
contours and velocity vectors on the meridional plane. The second and fourth 
rows show the field lines (lines of constant A v ) for the poloidal magnetic field at 
the same times as the first and third rows. The density contours are drawn 
for po/pmax(0) = io-°- 36i; -° 09 (j = 0-10). The field lines arc drawn for 



Oi/20 (t 



are the maximum and minimum values of A v respectively at the given time. 
Note that the field lines and the density contours show little change for t ^ 35P C , 
indicating that the star has settled down to an equilibrium state. 
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Figure 7. Evolution of (|_B J 9jf!|) (normalized to unity at t = 0). Note that this 
quantity drops toward zero at late time, indicating that the star is driven to a 
differentially rotating equilibrium state in which Q is constant along the magnetic 
field lines. 



constant along magnetic field lines. This is demonstrated in Fig. [7| which shows that 
the density-weighted average of the quantity B^djil approaches zero at late times. 
Since the rotation profile is adapted to the magnetic field structure, a stationary final 
state is reached which allows differential rotation and a nonzero toroidal field. 

44. Star B 2 

Both stars Bl and B2 are nonhypermassive. However, star Bl is ultraspinning, 
whereas B2 is normal. We evolve B2 with a seed magnetic field strength C = 
2.5 x 1CP 3 . Our simulation shows that this star evolves to a uniformly rotating 
configuration with little structural change. 

Figurc[S]shows the density contours and poloidal magnetic field at the initial time 
(t = 0) and at t = 46.5P C w 5.8^- We see that the density profile of the star does not 
change appreciably. This is not surprising since the main effect of the MHD processes is 
to redistribute the angular momentum inside the star. However, the rotational kinetic 
energy of star B2 is not very large (the initial = 0.040). Hence, the change of 

the centrifugal force inside the star as a result of angular momentum transport does 
not disturb the initial equilibrium significantly, unlike the cases of stars A and Bl. 

5. Discussion and Conclusions 

We have discussed in detail the evolution of magnetized HMNSs as first reported 
in [181 119) . In addition, we have performed simulations of two differentially rotating, 
but nonhypermassive, neutron star models with the same initial magnetic field 
geometry. These simulations have revealed a rich variety of behavior with possible 
implications for astrophysically interesting systems such as binary neutron star 
remnants, nascent neutron stars, and GRBs, where magnetic fields and strong gravity 
both play important roles. 
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Figure 8. Snapshots of the rest-mass density contours and poloidal magnetic field 
lines for star B2 at times t = and t = 46.5P C . The first row shows snapshots of 
the rest-mass density contours on the meridional plane. The second row shows the 
corresponding field lines for the poloidal magnetic field at the same times. The 
density contours are drawn for po/Pmax(0) = lo _t) - 36l— 009 [i = 0—10), where 
Pmax(0) is the maximum rest- mass density at t = 0. The field lines are drawn for 
A-ip = A^ m [ n + (A^ 
and 

meridional components of the velocity (which are zero initially) at t = 46.5P C are 
very small and so are not shown here. 



The hypermassive model considered here, star A, collapses to a BH due to the 
influence of the initial seed magnetic field. The early phase of evolution is dominated 
by magnetic winding, which proceeds until the back-reaction on the fluid becomes 
strong enough that the growth of the toroidal field ceases. This happens after roughly 
one Alfven time. After several rotation periods, we also see the effects of the MRI, 
which leads to turbulence. Eventually, the inner core of star A becomes unstable and 
collapses to a BH. Surrounding the BH, a significant amount of material remains in 
a magnetized torus which has been heated considerably by shocks resulting from the 
turbulent motions of the fluid. This final state consisting of a BH surrounded by a 
massive, hot accretion disk may be capable of producing highly relativistic outflows 
and a fireball (either through v — v annihilation or MHD processes) and is hence a 
promising candidate for the central engine of short-hard GRBs. This model predicts 
that such GRBs should accompany a burst of gravitational radiation and neutrino 
emission from the HMNS delayed collapse. 

The behavior of the nonhypermassive, ultraspinning star Bl under the influence 
of a seed magnetic field is quite different. Magnetic braking and the MRI operate in 
this model as well, leading to a mild contraction of the inner core and the expansion 
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of the outer layers into a high angular momentum torus-like structure. The final state 
consists of a fairly uniformly rotating core surrounded by a differentially rotating torus. 
The remaining differential rotation does not shear the magnetic field lines, so that the 
toroidal field settles down. On the other hand, the normal star B2 simply evolves to 
a uniformly rotating configuration. 

Two issues in particular warrant further study. The first is the scaling behavior 
of our solutions. We begin our simulations with a seed magnetic field which, though 
far too weak to be dynamically important, may be significantly larger than magnetic 
fields present in HMNSs formed through stellar collapse or a binary neutron star 
merger. We have demonstrated that, by varying the strength of the initial magnetic 
field through a factor of ~ 3 (See Fig. 0}, our evolution obeys the expected scaling 
during the magnetic winding phase, and the qualitative outcome of the simulations 
remains the same. However, since the MRI grows on a timescale ~ few x P c regardless 
of the initial magnetic field strength, it is possible that, for very weak initial fields, the 
effects of the MRI could dominate the evolution long before the effects of magnetic 
braking become important. In this case, the scaling of our numerical results with the 
Alfven time (relevant for magnetic winding) may break down. 

Another concern is the effect on our evolutions of relaxing the axisymmetry 
assumption. Rapidly and differentially rotating neutron stars may be subject to 
bar and/or one-armed spiral mode instabilities which could affect the dynamics 
(though star A was shown in Q to be stable against such instabilities, at least 
on dynamical timescales). Additionally, the development of the MRI in 2D differs 
from the 3D case [321 ■ Turbulence arises and persists more readily in 3D due to the 
lack of symmetry j^Sl OH] ■ However, McKinney and Gammie j^S] have performed 
axisymmetric simulations of magnetized tori accreting onto Kerr BHs and have found 
good quantitative agreement with the 3D results of De Villiers and Hawley on 
the global features of the evolution. Though simulations in full 3D will eventually be 
necessary to capture the full behavior of magnetized HMNSs, the 2D results presented 
here likely provide (at least) a good qualitative picture. 
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